# -*- coding: utf-8 -*-
"""
Created on Thu Apr 19 09:28:07 2012

@author: Daniele
"""

import numpy as np
import math as mat
import random as rd
import matplotlib.pyplot as plt
from pylab import cm
from numpy import exp



def Funcaov(ent,p):
# p4 = x =ent[p][0]
# Q24 = y  = ent[p][1]
# Q34 = z  = ent[p][2]

    f= np.zeros([D],float) #criando o vetor f
    p2=156.6-0.00752*ent[p][1]**2 #
    p3=117.1-0.00427*ent[p][2]**2 #
    D34=2.067/12.0 #
    D24=1.278/12.0 #
    D45=2.469/12.0 #
    fF=0.015 #
    pi=3.1416 #    
    k45=2*fF*145.0/((60.0*7.48)**2*(pi*D45**2.0/4.0)**2*D45) #
    k24=2*fF*125.0/((60.0*7.48)**2*(pi*D24**2.0/4.0)**2*D24) #
    k34=2*fF*125.0/((60.0*7.48)**2*(pi*D34**2.0/4.0)**2*D34) #    
    
    f[0]=70.0*32.3-ent[p][0]*(144*32.2/62.35)+k45*(ent[p][1]+ent[p][2])**2     
    f[1]=(ent[p][0]-p2)*(144*32.2/62.35)+k24*ent[p][1]**2
    f[2]=(ent[p][0]-p3)*(144*32.2/62.35)+k34*ent[p][2]**2
#    if (p==49):
#        print f[0]
#        print f[1]
#        print f[2]
    
    return (f)            



def minimo(x):
    
    minim= x[0]
    
    for i in range(1..len(x)):
        if(x(i)<minim):
            minim=x(i)    
    return(minim)

def maximo(x):
    
    maximo= x[0]
    
    for i in range(1..len(x)):
        if(x(i)>maximo):
            maximo=x(i)
    return(maximo)

def disteucli(ent):
    somat=0
    somat=(ent[0]*ent[0])+(ent[1]*ent[1])+(ent[2]*ent[2]);
    return (np.sqrt(somat))    


    
    
num_gen=200
num_vag=25#número de vagalumes..
D=3 #dimensão do problema
minimo1=0.0
maximo1=0.0
t=0

x= np.zeros([num_vag,D],float) #armazenando os dados no vetor
#print x 

int_lumx = np.zeros(num_vag,float)
r= np.zeros([num_vag,num_vag],float)
u= np.zeros([num_vag,num_vag],float)

minx = 56
maxx = 58

miny = 50
maxy  = 52

minz = 91
maxz  = 93

temp= np.zeros([D],float) #armazenando os dados no vetor

for i in range(num_vag):    
    x[i][0] = np.random.uniform(minx,maxx) #armazenando os dados no vetor
    x[i][1] = np.random.uniform(miny,maxy) #armazenando os dados no vetor
    x[i][2] = np.random.uniform(minz,maxz) #armazenando os dados no vetor
    temp=Funcaov(x,i)
    int_lumx[i]= disteucli(temp)

#print x
#print int_lumx
 
alfa=0.1
coef_ab=1
bzero=1
#intensidade luminosa....(aplicar os numero de vagalumes
# na função...

a1=0.0
b1=0.0
while (t<num_gen):
    t=t+1
    for i in range(num_vag):
        for j in range(num_vag):
            #if  (int_lumx[i]<int_lumx[j]): # busca de minimo
            if  (abs(int_lumx[i])>abs(int_lumx[j])): # busca de raiz 
                r[i][j]=np.sqrt ( (x[j,0]-x[i,0])**2 + (x[j,1]-x[i,1])**2 + (x[j,2]-x[i,2])**2)
                beta=bzero*np.exp(-1.0*(coef_ab* r[i][j]**2))
                x[i,0]= x[i,0]+beta*(x[j,0]-x[i,0])+(alfa*(np.random.uniform(0.0, 1.0,1)-0.5))
                x[i,1]= x[i,1]+beta*(x[j,1]-x[i,1])+(alfa*(np.random.uniform(0.0, 1.0,1)-0.5))
                x[i,2]= x[i,2]+beta*(x[j,2]-x[i,2])+(alfa*(np.random.uniform(0.0, 1.0,1)-0.5))
                temp=Funcaov(x,i)
                int_lumx[i]= disteucli(temp)


melhor = abs(int_lumx[0]);
for j in range(num_vag):
    if ( abs( int_lumx[j] ) < melhor ):
        melhor=abs(int_lumx[j])
        posmelhor=j
        
print x[posmelhor,0]
print x[posmelhor,1]
print x[posmelhor,2]
print Funcaov(x,posmelhor)


#print x 
#print y 
#print int_lumx

'''
pop = []
for i in range(num_vag):
    tmp = []
    tmp.append(x[i])
    tmp.append(y[i])
    pop.append(tmp[:])
    
    
def getCol(pop,col):
    z = []
    tam  = len(pop)
    for j in range(0,tam):
        z.append(pop[j][col])
    return z

def CurvaNivel(pop,l):
    plt.figure() # Create a new figure window
    xlist = np.linspace(-l, l, 100) # Create 1-D arrays for x,y dimensions
    ylist = np.linspace(-l, l, 100)
    X,Y = np.meshgrid(xlist, ylist) # Create 2-D grid xlist,ylist values
    Z = -(np.abs(X)+np.abs(Y))*np.exp(-(X**2 + Y**2)) # Compute function values on the grid
    x = getCol(pop,0) 
    y = getCol(pop,1)
    plt.contourf(X, Y, Z,cmap=plt.cm.jet)
    plt.colorbar() # adding the colobar on the right
    plt.scatter(x,y,marker='o',c='w',s=50)
    plt.axis([-l, l, -l, l]) # Change axis limits
    plt.title('$z=-(|x|+|y|) e^{-(x^2+y^2)}$')
    plt.xlabel('$x_1$')
    plt.ylabel('$x_2$')
    plt.show()

CurvaNivel(pop,D)


def z_func(x,y):
 return -(abs(x)+abs(y))*exp(-(x**2+y**2))

x = np.linspace(-5,5,100)
y = x
[X,Y] = np.meshgrid(x,y )
Z = z_func(X, Y)

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, 
                      cmap=cm.RdBu,linewidth=0, antialiased=False)

fig.colorbar(surf, shrink=0.5, aspect=5)
plt.title('$z=-(|x|+|y|) e^{-(x^2+y^2)}$')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.show()


'''
                           
#                if ( x[i] > x[j]): 
#                    a1=np.random.uniform(x[j],x[i],1)
#                else:
#                    a1=np.random.uniform(x[i],x[j],1)
                    

#                if ( y[i] > y[j]): 
#                    b1=np.random.uniform(y[j],y[i],1)
#                else:
#                    b1=np.random.uniform(y[i],y[j],1)
  